function [pointRoot] = bisectionMethod(range, tol, maxIter, typeTol)
%BISECTIONMETHOD Function that determines the root of a linear function
%using the Bisection Method.
%   This algorithm was adapted from Chapter 2.1 of Numerical Analysis (8th
%   edition, Burden & Faires)
%   
%   *Input Parameters*
%   range: [a b] where a<b, a and b are real numbers.
%   tol: Tolerance limit. 
%   maxIter: Maximum number of iterations
%   typeTol: String that defines how the tolerance condition is evaluated.
%            See tolCondition function in this file for available options.
%
%   *Default options* 
%   tol = 10^-5
%   maxIter = 100
%   typeTol = 'rangeTol' -> (pointB-pointA)/2
%
%   *Function to evaluate root is defined by myFunc at the bottom of 
%   this file.*
%
%   ***If no root is found after maxIteration, the midpoint for maxIter-1 
%   approximation is returned.***
%
%   ***********************************************************************
%   Author: Mathieu Boudreau, BSc, MSc, PhD Candidate (BME)
%   Institute: Montreal Neurological Institute, McGill University
%   Contact: mathieu.boudreau2 (at) mail.mcgill.ca
%   Date: July 15th 2014
%   ***********************************************************************


%% Set default input conditions, if required.
%
if nargin < 4
    typeTol = 'rangeTol';
end

if nargin < 3
    maxIter=100;
end

if nargin < 2
    tol=10^-5;
end


%% Organise input data
%
pointA=range(1);
pointB=range(2);

if sgn(myFunc(pointA))*sgn(myFunc(pointB))>=0
    error('Stop. Function values at both end of ranges are of the same sign, or one of them is equal to 0.')
end

%% Initialize initial conditions
%
iter=1;
funcA=myFunc(pointA);
oldPointMid = []; % Needed for 'absMidDiff' and 'relMidDiff' typeTol options.

%% Run Bisection Algorithm
%
while iter<maxIter
    % Set new midpoint, evaluate function at midpoit.
    pointMid = pointA+(pointB-pointA)/2;
    funcMid=myFunc(pointMid);
    
    if ~((strcmp(typeTol, 'absMidDiff') || strcmp(typeTol, 'relMidDiff')) && (iter == 1)) % These tolerance condition require at least two iteration
        
        % Verify if solution found, or if interval is smaller than tolerance
        if funcMid==0
            disp('Function evaluated at "exact" root point.')
            pointRoot = pointMid;
            break
        elseif tolCondition([pointA pointB], pointMid, oldPointMid, funcMid, typeTol)<tol
            disp('X interval is smaller than tolerance, point returned is middle of interval.')
            pointRoot = pointMid;
            break
        end
    end
    
    % Increase interval
    iter=iter+1;
    
    % Set new interval
    if sgn(funcA)*sgn(funcMid)>0
        pointA=pointMid;
        funcA=funcMid;
    else
        pointB=pointMid;
    end
    
    oldPointMid=pointMid;
end

if iter==maxIter
    disp('Function terminated after the maximum number iterations. No solution found within tolerance.')
    pointRoot = pointMid; 
end

end

function tolCalc=tolCondition(currentRange,currMidPoint,oldMidPoint,funcMid,typeTol)
    %% This function chooses and evaluates the tolerance condition. 
    % One of each flavour, take your pick! 
    %
    % (Burden and Faires suggest rangeTol as the best one for this problem, 
    % see discussion in Section 2.1 of their Numerical Analysis book (8th 
    % edition)
    
    switch typeTol
        case 'rangeTol'        
            tolCalc=(currentRange(2)-currentRange(1))/2;
        case 'absMidDiff'
            tolCalc=abs(currMidPoint-oldMidPoint);
        case 'relMidDiff'
            tolCalc=abs(currMidPoint-oldMidPoint)/abs(currMidPoint);
        case 'funcValue'
            tolCalc=abs(funcMid);
    end
    
end

function y = myFunc(x)
   %% This is the function that we are looking for roots in. Change it to match your problem!
   % 
   
   y=x-2^-x;
   
end